function [nx,a,sa,b,sb,r,sy] = wlinfit(xdat,ydat,eydat,ip)
%WLINFIT - least square straight line fit
%
%   [n,a,sa,b,sb,r,sy] = WLINFIT(xdat,ydat)
%   [n,a,sa,b,sb,r,sy] = WLINFIT(xdat,ydat,eydat)
%   [n,a,sa,b,sb,r,sy] = WLINFIT(xdat,ydat,eydat,IP)
%
% Performs the non-weighted or weighted least square straight line fit    
%
%        y =  a + b*x
%
% Input variables:
%   xdat : (mandatory input argument) x data, whose uncertainties are 
%          assumed to be negligible;
%   ydat : (mandatory input argument) y data;
%   eydat: (optional input) if eydat is undefined, empty, or is not a valid
%          vector of uncertainties (see below), a non-weighted data fitting 
%          to a straight line is carried out.
%          If eydat is defined and has the same size of ydat, a weighted 
%          fitting is carried out by using eydat as vector of uncertainties 
%          on y; 
%   IP   : if IP is undefined, empty, zero or false, no constraints about 
%          the intercept a are introduced.
%          If IP is true, the constraint a = 0 is introduced. 
%
% Output variables:
%   a, b  : straight line coefficients;
% 	sa, sb: standard deviation of a and b respectevely;
%	r     : linear correlation coefficient;
% 	sy    : standard deviation of mean(y).

% G. Teza, 2002, 2008, 2021 

xdat = xdat(:); % to have column vectors 
ydat = ydat(:);
if nargin < 3 || isempty(eydat)
    eydat = ones(size(ydat));
end

nx = numel(xdat); 
ny = numel(ydat);
if nx ~= ny
    error('The lengths of input vectors are different'); 
end
if nx <= 2
    error('Not enough input data; at least three experimental points must be available'); 
end 
if max(abs(xdat)) == 0
    error('x is the null vector'); 
end
if nargin < 2
    error('Not enough input variables'); 
end
if nargin < 4 || isempty(ip)
    ip = 0; 
end
ip = logical(ip);

I = find(eydat < eps*10^3);
if ~isempty(I)  % removal of points having invalid uncertainty
    xdat(I) = []; 
    ydat(I) = []; 
    eydat(I) = []; 
    nx = numel(xdat);
end
molty = ones(1,nx);  % for the dot product 

w = 1./(eydat.^2);  
sumw = molty*w;    
xm  = molty*(w.*xdat); 
x2m = molty*(w.*(xdat.^2));
ym  = molty*(w.*ydat); 
y2m = molty*(w.*(ydat.^2));
xym = molty*(w.*xdat.*ydat);
   
if ~ip
    denom  = sumw*x2m-xm^2; 
    denom2 = sumw*y2m-ym^2;
    numa   = x2m*ym-xm*xym;
    nullor = 1; 
    numb   = sumw*xym-xm*ym;
else
  	denom  = sumw*x2m; 
    denom2 = sumw*y2m;
  	numa   = 0; 
    nullor = 0; 
  	numb   = sumw*xym;
end   

a  = numa/denom; b = numb/denom;
yc = a+b*xdat;  

err  = ydat-yc; 
sumq = molty*(err.^2);
sy   = sqrt(sumq/(nx-2));
fatco = sy/sqrt(nx);

sa = fatco*nullor*sqrt(x2m/denom); 
sb = fatco*sqrt(1/denom);
if abs(denom2) > eps*10^3
    r = numb/sqrt(denom*denom2);
else
    r = 1;  % In such a case, an horizontal straight line is obtained
end